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: ABSTRACT 

cn 

The vertical structure of black hole accretion disks in which radiation domi- 

^ ^ nates the total pressure is investigated using a three-dimensional radiation-MHD 

Q> \ calculation. The domain is a small patch of disk centered 100 Schwarzschild radii 

^ I from a black hole of 10* Mq, and the stratified shearing-box approximation is 

CN \ used. Magneto-rotational instability converts gravitational energy to turbulent 

^ I magnetic and kinetic energy. The gas is heated by magnetic dissipation and by 

O I radiation damping of the turbulence, and cooled by diffusion and advection of ra- 

' diation through the vertical boundaries. The resulting structure differs in several 

^ ' fundamental ways from the standard Shakura-Sunyaev picture. The disk consists 

' of three layers. At the midplane, the density is large, and the magnetic pressure 



^ ■ and total accretion stress are less than the gas pressure. In lower-density surface 

layers that are optically thick, the magnetic pressure and stress are greater than 
the gas pressure but less than the radiation pressure. Horizontal density varia- 



;h ' tions in the surface layers exceed an order of magnitude. Magnetic fields in the 

regions of greatest stress are buoyant, and dissipate as they rise, so the heating 
rate declines more slowly with height than the stress. Much of the dissipation oc- 
curs at low column depth, and the interior is cooler and less radiation-dominated 
than in the Shakura-Sunyaev model with the same surface mass density and flux. 
The mean structure is convectively stable. 

Subject headings: accretion, accretion disks — instabilities — MHD — radiative 
transfer 
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1. INTRODUCTION 

Black holes in X-ray binary star systems and active galactic nuclei accrete material 
through surrounding disks of gas. When the accretion rate is greater than about 1% of the 
Eddington value, the radiation pressure in the central regions of the disk may exceed the 
gas pressure (Shakura & Sunyaev 1973). The disk structure and the radiation emerging 
are governed by removal of angular momentum from the gas, and conversion of released 
gravitational energy to photons. In the standard model, the angular momentum transfer 
and dissipation processes are unspecified. Both are assTimed to occur at rates proportional 
to the height-integrated total pressure, and the heating is balanced by diffusion of photons 
to the disk surfaces (Shakura & Sunyaev 1973). The model is unstable to perturbations in 
mass accretion rate (Lightman & Eardley 1974) and heating rate (Shakura & Sunyaev 1976), 
and to vertical convective displacements (Bisnovatyi-Kogan & Blinnikov 1977). Any one of 
these instabilities might substantially modify the disk structure. However, if the stresses are 
magnetic in nature, buoyancy of the fields may lead to stresses scaling with the gas pressure 
alone, rather than the total pressure (Sakimoto & Coroniti 1989). In this case the viscous 
and thermal instabihties are absent. The instabilities may also be eliminated if part of the 
accretion energy is dissipated outside the disk, in a hot corona (Svensson & Zdziarski 1994). 
Observed properties that are not understood using the standard Shakura- Sunyaev model 
include the X-ray emission (Elvis et al. 1978) and UV-optical spectra (Blaes 2004) of active 
galactic nuclei, and the steady disk luminosities of some stellar-mass black hole X-ray binary 
systems accreting at 1% to 50% of the Eddington rate (Gierlihski & Done 2004). 

The most likely physical mechanism of angular momentum transfer in radiation-domin- 
ated disks is magnetic stresses. The magneto-rotational instability or MRI (Balbus & Hawley 
1991) leads to turbulence in which the energy of differential orbital motion is converted to 
magnetic fields and gas motions (Hawley, Gammie, & Balbus 1996). In isothermal MHD 
simulations in which heating and cooling are assumed to balance, magnetic fields generated 
in the turbulence are partly expelled to form magnetized coronae above and below the disk 
(Miller & Stone 2000). 

Disk structure may be affected also by the locations and rates of heating and cool- 
ing. Candidate physical heating mechanisms include resistive dissipation of magnetic fields, 
microscopic viscous dissipation of gas motions, and radiative damping of compressible tur- 
bulence (Agol & Krohk 1998). Coohng processes may include radiation diffusion, convection 
(Agol et al. 2001), and photon bubble instability (Gammie 1998). Here I examine the ef- 
fects of heating and cooling on the vertical structure of radiation-dominated disks, using a 
radiation-MHD simulation. 
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2. DOMAIN AND METHODS 

The domain is a small patch of disk, centered 100 Schwarzschild radii Rs = 2GM/c^ 
out from a black hole of M = 10^ Mq, and co- rotating at the central Keplerian orbital 
frequency VLq. The local shearing-box approximation is used, and Cartesian coordinates 
(x, I/, z) correspond to distance from the origin along the radial, orbital, and vertical direc- 
tions, respectively (Hawley, Gammie, & Balbus 1995). The domain extends l.bRs along the 
radial direction, Q)Rs along the orbit, and QRs either side of the midplane, and is divided 
into 32 X 64 X 256 zones. 

The equations solved and numerical method differ from Turner et al. (2003) in including 
the vertical component of the gravity of the black hole by a term —pVt^z in the equation of 
motion. No viscosity of the Shakura-Sunyaev type is used. The frequency-averaged equations 

of radiation MHD (Mihalas & Mihalas 1984; Stone, Mihalas, & Norman 1992) arc integrated 
in the fiux-limitcd diffusion (FLD) approximation (Levermore & Pomraning 1981), using 
the Zeus code (Stone & Norman 1992a,b) with its FLD module (Turner & Stone 2001). 
Opacities due to electron scattering and free-free processes are included. The equations are 
closed with a 7 = 5/3 ideal-gas equation of state. 

Several dissipation mechanisms may act. A total energy scheme is used during the 
magnetic field update, so that numerical losses of magnetic field are captured as gas heat. 
The gas internal energy may also increase through shock compression and associated artificial 
viscosity. The radiation energy increases when photons diffuse from compressed regions, 
irreversibly extracting part of the work done in compression (Agol & Krolik 1998; Turner, 
Stone, & Sano 2002). Gas and radiation remain near mutual thermal equilibrium while 
heating, owing to free-free absorption and emission. For a steady state, the total dissipation 
must be balanced by diffusion of radiation, and advection of gas, radiation, and magnetic 
energy through the vertical boundaries. 

The azimuthal boundaries are periodic, the radial boundaries shearing-periodic, and 
the vertical boundaries allow outflow but no inflow. The radiation flux into each vertical 
boundary zone is set equal to that between the two adjacent active zones, except that 
the fiux is fixed at zero if needed to prevent radiation energy from entering the domain. 
Updated boundary fiuxes are estimated using photon diffusion coefficients from the previous 
timestep. A density floor is applied throughout. In zones where the density falls below 0.2% 
of the initial midplane value, mass is added to bring the density up to the floor level. The 
corresponding minimum optical depth per zone in the vertical direction is 24. 
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3. INITIAL STATE 

The surface mass density in the calculation sets the total optical depth, and varies from 
its initial value only through outflows and the density floor. The net vertical magnetic flux 
affects the accretion stress (Hawley, Gammie, & Balbus f995), and is time-constant in the 
shearing-box approximation. Other aspects of the initial condition likely have little effect on 
the outcome, since the structure is to be determined by accretion stresses, dissipation, and 
cooling. 

The initial condition is a Shakura & Sunyaev (1973) model accreting at 10% of the 
Eddington value for luminous efficiency 0.1. In constructing the initial state only, the ratio 
a of stress to total pressure is set to 0.01. The resulting surface density is 1.1 x 10^ g cm~^. 
The half-thickness B. of the Shakura-Sunyaev model is \Rs- The domain outside the model 
is filled with an ambient medium of the floor density. Since gravity increases with height, 
while radiation flux is independent of height outside the Shakura-Sunyaev model, the ambient 
medium is out of hydrostatic balance and falls towards the midplane when the calculation 
starts. 

The magnetic fleld is given zero net vertical flux, so that outflows can completely de- 
magnetize the domain. The starting conflguration is an azimuthal flux tube of circular 
cross-section, with radius 0.75H. The axis is offset from domain center by +0.1H in x and 
+0.1H in z. Field strength in the tube is uniform at 2660 Gauss, corresponding to 4% 
of initial midplane gas plus radiation pressure. The tube is twisted about its axis, giving 
maximum poloidal component 661 Gauss. The maximum vertical MRI wavelength of 8 grid 
zones is adequately resolved. The domain-mean radial field is zero, and the mean azimuthal 
field of 159 Gauss is much less than values that develop later. The calculation is begun 
with a small random poloidal velocity in each grid zone. The maximum amplitude of each 
velocity component is 1% of the midplane radiation sound speed. 



4. RESULTS 

During the first two orbital periods, the fiux tube is stretched radially by MRI, and its 
upper parts are lifted to the surface of the Shakura-Sunyaev model by magnetic buoyancy. 
By four orbits, the tube is torn apart by MRI and spread throughout the region within Rg 
of the midplane. Gravitational energy is released faster than assumed in the initial Shakura- 
Sunyaev model, and the magnetized region expands into the ambient medium. After 13 
orbits, magnetic fields are found throughout the domain. 
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4.1. Three Layers 

The horizontally- averaged structure present from 13 orbits on may be divided into three 
layers (Figure 1). In a dense layer within about 1.5Rs of the midplane, mean magnetic 
pressure and total accretion stress are less than gas pressure. In lower-density surface layers, 
magnetic pressure and accretion stress are greater than gas pressure, but less than the sum 
of gas and radiation pressures. The average structure is marginally convectively stable in the 
midplane layer. In the surface layers, the gas and radiation Brunt- Vaisala frequencies are 
about equal to the orbital frequency Qq, and are real, indicating hydrodynamic convective 
stability. 

Magnetic energy is produced fastest near the heights of greatest accretion stress. The 
fields are buoyant and rise towards the boundaries (Figure 2). The rise speed is approximately 
the Alfven speed, as expected for magnetized regions exchanging heat rapidly with their 
surroundings (Parker 1975). The fields dissipate numerically while rising, and the mean fiux 
of magnetic energy through the boundaries between 15 and 45 orbits is 6% of the mean 
radiative fiux. The horizontally-averaged magnetic pressure is dominated by the azimuthal 
component throughout, and in most cases the sign of the azimuthal field alternates between 
consecutive rising regions. The regular pattern seen in figure 2, repeating about every 7 
orbits, may result in part from the limited domain size. Rising regions in most cases fill the 
horizontal extent of the box. In real disks, adjacent buoyant regions may interact. 

The photons may partly decouple from the turbulence, in radiation-supported disks 
accreting via MRI, as the distance radiation diffuses per orbit is about equal to the vertical 
MRI wavelength (Turner et al. 2003). Large density variations accompany loss of radiation 
pressure support if magnetic pressure exceeds gas pressure (Turner, Stone, & Sano 2002). 
The results of the present calculation are consistent with this picture. Density variations on 
horizontal x — y planes exceed an order of magnitude near the upper and lower boundaries, 
where the radiation diffusion distance is half the RMS vertical MRI wavelength, and the 
mean magnetic pressure is more than 30 times the gas pressure. By contrast, horizontal 
density variations within Rs of the midplane are only a factor two. The average diffusion 
distance is equal to the RMS vertical MRI wavelength, but the mean magnetic pressure is 
less than the gas pressure, so magnetic forces squeeze the gas less. At the height 3.3Rs 
where vertical magnetic pressure peaks, the diffusion length is 0.3 times the vertical MRI 
scale. Although radiation is better-coupled to gas motions, the magnetic pressure is greater 
than the gas pressure. Density variations are intermediate between those at midplane and 
boundaries. 

The vertical variation in magnetic field strength may be related also to the diffusion 
of photons through the turbulence. Linear MRI grows slowly if magnetic pressure exceeds 
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Fig. 1. — Simulation results averaged horizontally and over time from 15 to 45 orbits. The 
three layers may be seen at top. Curves show density (solid, left scale) and magnetic pressure 
(dashed, right scale) versus height. The surface layers are optically thick. The optical depth 
between each magnetic pressure peak and the nearby boundary is 1.5 x 10^. At bottom are 
plotted the pressures due to radiation (solid), gas (dotted) and magnetic fields (dashed), and 
the total accretion stress (dot-dashed). 
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the pressure resisting compression (Blaes & Balbus 1994; Blaes & Socrates 2001). Near 
the midplane in the simulation, where fluctuations are due largely to the competition of 
magnetic with gas pressure forces, the turbulence saturates with magnetic pressure less than 
gas pressure. Around 3.3Rs, where squeezing is resisted by both radiation and gas pressures, 
the magnetic pressure is greater than the gas pressure. Near the vertical boundaries, the 
rate of arrival of magnetic fields from the interior exceeds the local rate of field generation. 

The vertical magnetic pressure profile is similar to those in isothermal MHD simulations 
of gas-supported disks by Miller & Stone (2000). However the magnetized surface layers 
here arc optically thick, and lie inside the disk photosphere. Horizontally-averaged magnetic 
pressure is everywhere less than gas plus radiation pressure. The total accretion stress at the 
midplane is half the maximum in the surface layers, whereas in the Miller & Stone (2000) 
fiducial calculation, total stress in the gas-dominated region varies irregularly with height. 

4.2. Comparison with Shakura-Sunyaev Model 

The averaged simulation results may be compared against the Shakura-Sunyaev model 
with the same surface density and radiation fiux. The model has a = 0.0013 and accretion 
rate 64% of the Eddington value. Density is almost independent of height in the model 
interior, while in the simulation the density is centrally-concentrated (Figure 3). In the 
Shakura-Sunyaev picture, it is often assumed that dissipation is locally proportional to den- 
sity. By contrast, much of the dissipation in the simulation occurs in regions of low density. 
Photons escape easily from the surface layers, so the interior is cooler than in the Shakura- 
Sunyaev model. The midplane ratio of radiation to gas pressure is 160 in the model, and 
14 in the simulation. The thermal time in the Shakura-Sunyaev model is 370 orbits, while 
in the simulation the cooling time, or ratio of total energy content to cooling rate, is just 
19 orbits. Dissipation falls off more slowly with height than the stress, as fields tend to rise 
between generation and dissipation. 

4.3. Thermal Balance 

The radiation content of the fiow is maintained through accretion. Between 15 and 45 
orbits, the energy released is 157% of the mean domain-integrated total energy, while the 
total energy is almost constant, increasing just 2%. The source of power, the differential 
rotation, is transformed into turbulent magnetic and kinetic energy by magnetic stresses. 
The turbulence is converted to heat mainly through numerical losses of magnetic fields. 
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Radiation damping contributes 29% of the heating within Rs of the midplane, where it is 
readily distinguished from expansion associated with magnetic buoyancy. Artificial viscous 
heating, mostly in shocks, makes up 12% of the overall dissipation. The dissipation is 
approximately balanced by cooling. The main cooling processes arc diffusion and advection 
of radiation through the vertical boundaries, with diffusion carrying two-thirds of the energy 
flux. Mass lost through the boundaries is balanced by mass added to maintain densities 
above the floor, and total mass increases by less than 1% from 15 to 45 orbits. The energy 
flux is 73% of the energy input, and the remainder disappears partly through numerical 
losses of kinetic energy. A complete energy-conserving scheme may be useful for future 
calculations. 

Starting at 50 orbits, the disk heats, then cools. Total energy increases roughly linearly 
by a factor 2.5 up to 90 orbits, decreases linearly until 145 orbits, then is again steady. The 
run ends at 170 orbits, or 8 simulation thermal times after flrst saturation. The heating rate 
varies, while the diffusion cooling rate is approximately time-constant. No exponentially- 
growing thermal instability is seen. However, during the hot period, the disk expands, and 
losses through the vertical boundaries reduce the mass by almost half. The absence of 
runaway heating may be due to the mass loss, rather than internal thermal stability. 

5. CONCLUSIONS 

A patch of radiation-dominated accretion disk lying 100 Schwarzschild radii from a 10^ 
M0 black hole is simulated including physical processes of energy release, dissipation, and 
cooling. The patch has surface density 10^ g cm^^, and its net vertical magnetic flux is 
zero. The structure that develops has density greatest at the midplane, stresses greatest in 
optically-thick surface layers, and dissipation found throughout but intermittent in time and 
space. The stresses result from magnetic forces. The flelds in the surface layers are buoyant. 
Heating occurs mainly by numerical dissipation of the flelds, and physical radiation damping 
of the turbulence, while cooling occurs by diffusion and advection of radiation through the 
boundaries. The vertical structure may prove to vary with mass of the central body, location 
in the disk, surface density, and magnetic flux. 

The time-averaged structure in the simulation is hydrodynamically convectively stable. 
Thermal instability may be prevented by outflow through the boundaries. Viscous stability 
is not tested, as there is no net radial flow in the shearing box. Photon bubble instability 
may be present. The fastest linear modes (Blaes & Socrates 2003) in the time-averaged 
structure are found in the surface layers, have wavelengths shorter than the vertical grid 
spacing, and grow at O.9f2o- Resolved modes, with wavelengths greater than six grid zones. 
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are expected to grow more slowly than the MRI. No clear signs of these modes are found. 
The hkely direct effect of photon bubbles is to further reduce the coohng time. 

The methods used here were developed with Jim Stone and Takayoshi Sano. I benefited 
also from discussions with Julian Krolik and Omer Blaes. The work was supported by the 
U. S. National Science Foundation under grant AST-0307657. 
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Fig. 2. — Horizontally-averaged magnetic pressure versus height and time. The grey scale is 
linear from zero (white) to 3 x 10^ dyn cm~^ (black). 
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Fig. 3. — Simulation results (solid) compared with the Shakura-Sunyaev model (dotted) 
having the same surface density and radiation flux. Density is shown at top, total dissipation 
below. The simulation results are averaged horizontally, and over time from 15 to 45 orbits. 



